Solid helium at high pressure: A path-integral Monte Carlo simulation 
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Solid helium (^He and ''He) in the hep and fee phases has been studied by path-integral Monte 
Carlo. Simulations were carried out in the isothermal-isobaric (NPT) ensemble at pressures up to 
52 GPa. This allows one to study the temperature and pressure dependences of isotopic effects on 
the crystal volume and vibrational energy in a wide parameter range. The obtained equation of 
state at room temperature agrees with available experimental data. The kinetic energy, Ek, of solid 
helium is found to be larger than the vibrational potential energy, Ep. The ratio Ek/Ep amounts to 
about 1.4 at low pressures and decreases as the applied pressure is raised, converging to 1, as in a 
harmonic solid. Results of these simulations have been compared with those yielded by previous path 
integral simulations in the NVT ensemble. The validity range of earlier approximations is discussed. 



PACS numbers: 67.80.-s, 62.50.-fp, 65.40.De, 05.10.Ln 



I. INTRODUCTION 

Structural and thermodynamic properties of solid he- 
lium have been of continuous interest in condensed mat- 
ter physics because of its quantum nature and electronic 
simplicity. In fact, solid helium is in many respects an 
archetypal "quantum solid," where zero-point energy and 
associated anharmonic effects are appreciably larger than 
in most known solids. This gives rise to peculiar prop- 
erties, whose understanding has presented a challence 
for elaborated theories and modelling from a microscopic 
standpoint 0. Among these properties, the behavior of 
condensed helium at high density has received much at- 
tention. In fact, diamond-anvil-cell and shock-wave ex- 
periments have allowed to study the equation of state 
(EOS) of solid helium up to pressures on the order of 
50 GPa 0,0, 0- In last years, the effect of pressure on 
heavier rare-gas solids has also been of interest for both 
experimentalists 0, |^ and theorists 0, H, IS ^ • 

Anharmonic effects in solids, and in solid helium in 
particular, have been traditionally studied by using theo- 
retical techniques such as quasiharmonic approximations 
and self-consistent phonon theories 11]. In more recent 
years, the Feynman path-integral formulation of statisti- 
cal mechanics been exploited to study ther- 
modynamic properties of solids at temperatures lower 
than their Debye temperature Qd, where the quantum 
character of the atomic nuclei is relevant. Monte Carlo 
sampling applied to evaluate finite-temperature path in- 
tegrals allows one to carry out quantitative and nonper- 
turbative studies of highly-anharmonic effects in solids 

Q 

The path-integral Monte Carlo (PIMC) technique has 
been used to s tudy several properties of solid helium 
[1 13 111 IB llTl. as well as heavier rare-gas solids 
[isL lig. 1m 12111221 . For helium, in particular, this 
method has predicted kinetic-energy values 0| and 
Debye- Waller factors |2^ in good agreement with data 
derived from experiments |3 US- PIMC simulations 



were also employed to study the isotopic shift in the he- 
lium melting pressure 0, The EOS of solid he- 
lium at T = has been studied by diffusion Monte Carlo 
in a wide density range (down to a molar volume of 15 
cm'^/mol), using an accurate interatomic potential |26| . 
This has been done by Chang and Boninsegni sA, fi- 
nite temperatures for both solid and liquid helium, by 
using PIMC simulations with several interatomic poten- 
tials, and for molar volumes down to 2.6 cm^/mol. These 
authors suggested that the use of effective potentials in- 
cluding two- and three-body terms alone can be insuffi- 
cient to reproduce the EOS of condensed helium, in the 
pressure range experimentally accessible at present. 

In this paper, we study the effect of pressure on solid 
■^He and ^He by PIMC simulations. We employ the 
isothermal-isobaric (NPT) ensemble, which allows us to 
consider properties of these solids along well-defined iso- 
bars. The interatomic interaction is described by a com- 
bination of two- and three-body terms, which are directly 
included in the simulations. This permits us to check 
results of earlier simulations, where the effect of three- 
body terms was taken into account in a perturbative way 
|ig.ll7l |. By comparing results for '^He and ^He, we anal- 
yse isotopic effects on the vibrational energy and crystal 
volume of solid helium. 

The paper is organized as follows. In Sec. II, the com- 
putational method is described. In Sec. Ill we present 
results for the equation of state, vibrational energy, and 
isotopic effects on the lattice parameters. Finally, Sec. IV 
includes a discussion of the results and the conclusions. 



II. METHOD 

Equilibrium properties of solid helium in the face- 
centred cubic (fee) and hexagonal close-packed (hep) 
phases have been calculated by PIMC simulations in the 
NPT ensemble. Most of our simulations were performed 
on supercells of the fee and hep unit cells, including 500 
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and 432 helium atoms respectively. To check the conver- 
gence of our results with system size, some simulations 
were carried out for other supercell sizes, and it was found 
that finite-size effects for N > 400 atoms are negligible 
for the quantities studied here (they are smaller than the 
error bars). 

Helium atoms have been treated as quantum parti- 
cles interacting through an effective interatomic poten- 
tial, composed of a two-body and a three-body part. For 
the two-body interaction, we have taken the potential de- 
veloped by Aziz et al. 23} (the so-called HFD-B3-FCI1 
potential). For the three-body part we have employed 
a Bruch-McGee-type potential p^ . [2^ . which includes 
the exchange three-body interaction and a triplc-dipole 
Axilrod- Teller interaction. For most of the simulations 
presented below, the parameters employed for the three- 
body terms were those given by Loubeyre but with 
the parameter A in the attractive exchange interaction 
rescaled by a factor 2/3 (as suggested in Ref. J^J, and 
giving A = 20.43 au). 

In the path-integral formulation of statistical mechan- 
ics, the partition function is evaluated through a dis- 
cretization of the density matrix along cyclic paths, 
composed of a finite number L (Trotter number) of 
"imaginary-time" steps [T^ .13 ] . In the numerical sim- 
ulations, this discretization gives rise to the appearance 
of L replicas for each quantum particle. Thus, the prac- 
tical implementation of this method relies on an isomor- 
phism between the quantum system and a classical one, 
obtained by replacing each quantum particle by a cyclic 
chain of L classical particles, connected by harmonic 
springs with a temperature-dependent constant. Details 
on this computational method can be found elsewhere 

QIMIlil. 

Our simulations were based on the so-called "primi- 
tive" form of PIMC [s^. l3^ . For interatomic potentials 
including only two-body terms, effective forms for the 
density matrix have been developed, which allow one to 
reduce efficiently the Trotter number, thereby simplify- 
ing appreciably the calculation |l6j| . Such a simplification 
is not possible here, since we consider explicitly three- 
body terms in the simulations. Quantum exchange ef- 
fects between atomic nuclei were not considered, as they 
are negligible for solid helium at the pressures and tem- 
peratures studied here. (This should be valid as long as 
there are no vacancies and T is greater than the exchange 
frequency ~ 10~^ K 0.) The dynamic effect of the in- 
teractions between nearest and next-nearest neighbours 
is explicitly considered. The effect of interactions be- 
yond next-nearest neighbours is taken into account by a 
static-lattice approximation which was employed 

earlier in PIMC simulations of rare-gas solids j^, |2^ . We 
have checked that including dynamical correlations be- 
tween more distant atom shells does not change the re- 
sults presented below. For the energy we have used the 
"crude" estimator, as defined in Refs. H^,!!!]- All calcu- 
lations have been performed using our PIMC code [35l |. 
that has been employed earlier to study various types of 



solids ■25,131131133. 

Sampling of the configuration space has been carried 
out by the Metropohs method at pressures P < 52 GPa, 
and temperatures between 25 K and the melting temper- 
ature of the solid at each considered pressure. A simu- 
lation run proceeds via successive Monte Carlo steps. In 
each step, the replica coordinates are updated according 
to three different kinds of sampling schemes: (1) Sequen- 
tial trial moves of the individual replica coordinates; (2) 
trial moves of the center of gravity of the cyclic paths, 
keeping unaltered the shape of each path, and (3) trial 
changes on the logarithm of the volume of the simulation 
cell. For given temperature and pressure, a typical run 
consisted of 10'' Monte Carlo steps for system equilibra- 
tion, followed by 10^ steps for the calculation of ensemble 
average properties. Other technical details are the same 
as those used in Refs. p^lsl). 

To have a nearly constant precision for the simulation 
results at different temperatures, we have taken a Trot- 
ter number that scales as the inverse temperature. At a 
given T, the value of L required to obtain convergence 
of the results depends on the Debye temperature O^i of 
the considered solid (higher Qd needs larger L). Since 
vibrational frequencies (and the associated Debye tem- 
perature) increase as the applied pressure is raised, L has 
to be increased accordingly. Thus, we have taken LT = 
4000 K for ^He and 3000 K for "^He, which were found to 
be sufficient for simulations of the corresponding solids 
at the pressures considered here (P < 52 GPa). This 
means that, for solid "^He at T = 25 K, we have L = 160 
and then the computational time required to carry out 
a PIMC simulation for N = 500 helium atoms is equiva- 
lent to a classical Monte Carlo simulation of LN ~ 80000 
atoms (assuming the same number of simulation steps) . 



III. RESULTS 

A. Crystal volume 

Shown in Fig. 1 is the pressure dependence of the 
crystal volume for solid "^He at 300 K. Open symbols 
represent results of PIMC simulations with different in- 
teratomic potentials: (a) only two-body interactions with 
an Aziz- type potential |23 (squares); (b) two-body inter- 
actions as in Ref. |23| plus three-body terms as in Ref. 
(29i | (diamonds); (c) the same two-body potential and 
three-body interaction with the exchange part rescaled 
by 2/3, as proposed in Ref. 0| (circles). For compar- 
ison, we also present results derived from PIMC simu- 
lations in the NVT ensemble with the exchange 
interaction rescaled by the same factor 2/3 (triangles). 
Filled symbols indicate experimental results obtained by 
Mao et al. Q (filled circles) and Loubeyre et al. 
(filled squares). Our simulation data show that the in- 
teratomic potential (c) gives results for the equation of 
state of solid ''He close to the experimental data. The 
only consideration of two-body terms predicts, for a given 
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FIG. 1: Equation of state pressure-volume of hep ''He at 300 
K, as derived from PIMC simulations for different interatomic 
potentials, (a) Open squares: only two-body interactions with 
an Aziz- type potential |2^ : (b) open diamonds: two-body 
terms as in Ref. |27t and three-body interactions, as defined in 
Ref. [2^; (c) open circles: same potential as (b), but with the 
exchange three-body interaction rescaled by 2/3. Error bars 
of the simulation results are less than the symbol size. Dotted 
lines are guides to the eye. Open triangles show earlier results 
obtained from PIMC simulations in the NVT ensemble [T^ . 
with the attractive exchange interaction rescaled as in (c). 
Black symbols show experimental data obtained by Mao et 
al. Q (filled circles) and Loubeyre et al. 4] (filled squares). 



applied pressure, a crystal volume larger than the exper- 
imental one. On the contrary, consideration of both two- 
and three-body terms derived from ab initio calculations 
underestimates the volume of solid helium. This is in 
line with results obtained earlier from PIMC simulations 
in the NVT ensemble in Refs. These authors 

introduced the three-body terms not directly in the sim- 
ulations, but calculated in a perturbative way their con- 
tribution to thermodynamic averages from configurations 
obtained in PIMC simulations. In Fig. 1 (open circles 
and triangles) we observe that our NPT simulations give 
the same results as those obtained earlier in the NVT en- 
semble (17i| for pressures lower than 30 GPa. However, 
at higher pressures both sets of results differ appreciably, 
and for a given P the procedure employed in Ref. [l7j 
yields a volume larger than that obtained here. 

Our results with the effective interatomic potential (c) 
follow closely the experimental ones, even at the highest 
pressures considered here, although they seem to become 
lower than the later as pressure rises. However, taking 
into account the dispersion of experimental points and 
the error bars of the simulation data, differences between 
both sets of data are not enough to invalidate the accu- 
racy of the effective potential (c) in the pressure range 
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FIG. 2: Temperature dependence of the lattice parameter a 
of hep helium, derived from PIMC simulations at two different 
pressures: (a) 12 GPa and (b) 52 GPa. Squares and circles 
indicate results for *He and ^He, respectively. Error bars are 
smaller than the symbol size. Dashed lines are guides to the 
eye. 



studied here. This contrasts with the conclusions pre- 
sented by Chang and Boninsegni ^"^o argued that 
two- and three-body terms alone may be insufficent to 
reproduce quantitatively the EOS of condensed helium 
at pressures on the order of 50 GPa. In fact, for molar 
volumes smaller than 3 cm'^/mol, these authors found 
pressures larger than the experimental ones (see Fig. 1). 
In view of these results, in the remainder of the paper we 
will employ only the interaction potential (c) (with the 
exchange three-body part of Ref. '23| rescaled by 2/3). 
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FIG. 3: Isotopic effect on the lattice parameter a of hep 
helium, as obtained from PIMC simulations. Shown is the 
parameter x = (is ~ a4)/ffl4 as a function of temperature for 
three different pressures: squares, 12 GPa; circles, 32 GPa; 
triangles, 52 GPa. Error bars are on the order of the symbol 
size. Lines are guides to the eye. 



In Fig. 2 we present the temperature dependence of 
the lattice parameter a for hep helium at two pressures: 
(a) 12 GPa and (b) 52 GPa. Squares and circles repre- 
sent results of our PIMC simulations for ^He and ^He, 
respectively. Note the different vertical and horizontal 
scales in (a) and (b). As expected, the difference — 
between lattice parameters of '^He and *He at a given 
pressure decreases as the temperature is raised (at high 
T the solid becomes "more classical"). 

To quantify the isotopic effect on the linear dimen- 
sions of solid helium, we employ the parameter x — 
(as — a4)/a4, which measures the relative difference be- 
tween lattice parameters of solid ^He and ^He. This pa- 
rameter X is displayed in Fig. 3 as a function of tem- 
perature for three different pressures: 12 (squares), 32 
(circles), and 52 GPa (triangles). For each pressure, re- 
sults are shown for temperatures at which the considered 
solids were stable along the PIMC simulations. At low 
T, the parameter x is larger for lower pressure, and val- 
ues obtained for different pressures approach each other 
as the temperature rises. At T > 300 K, we find that 
X is larger for P = 52 GPa than for 32 GPa. This does 
not mean that the difference 03 — 04 is larger at 52 GPa, 
but is a consequence of the normalization by 04, which 
is clearly smaller at higher pressure. 

The difference 03 — 04 is largest at small pressures and 
low temperatures, where quantum effects are most promi- 
nent. For a given solid, quantum effects on the crystal 
size can be measured by the difference Aa = a — ad be- 
tween the actual lattice parameter a and that obtained 



for a "classical" crystal of point particles, ad. This dif- 
ference decreases for increasing atomic mass and temper- 
ature [Ullil. From our PIMC simulations at T = 25 
K and a relatively low pressure of 0.3 GPa, we found an 
increase in the linear size of solid "^He and ^ He of 8.2 and 
7.2% with respect to the classical crystal at zero tem- 
perature. For comparison, we note that the "zero-point 
expansion" for heavier rare gases causes a relative in- 
crease of 4.1% and 1.2% in the lattice parameter of fee 
Ne and Ar, respectively [38l| . 

B. Energy 

For a given interatomic potential, the internal energy 
of a solid, i?(y, T), at volume V and temperature T can 
be written as: 

E{Y, T) = E^iniV) + E^ii^iV, T) , (1) 

where E^i^{V) is the potential energy for the (classical) 
crystal at T = with point atoms on their lattice sites, 
and i?vib(^, T) is the vibrational energy. Since we are 
working here in the isothermal-isobaric ensemble, the vol- 
ume is implicitly given by the applied pressure. At finite 
temperatures, V changes with T due to thermal expan- 
sion, and for real (quantum) solids, the crystal volume 
depends on quantum effects, which also contribute to ex- 
pand the crystal with respect to the classical expectancy 
(mainly at low T). Thus, the volume V in Eq. (^) is an 
implicit function of P and T, i.e., V = V{P,T). 

The vibrational energy, Eyi\,{y,T), depends explicitly 
on both, V and T, and can be obtained by subtracting 
the energy i?inin(^) from the internal energy. In this 
way, path-integral Monte Carlo simulations allow one to 
obtain separately the kinetic, -Efc, and potential energy, 
iJp, associated to the lattice vibrations: ii^vib = Ek + Ep. 

In Fig. 4 we show the kinetic energy of solid helium 
(hep phase) as a function of temperature for P = 12 
and 52 GPa. Open symbols correspond to "'He and filled 
symbols to ■^He. At 52 GPa, the low-temperature kinetic 
energy of solid ''He is found to be 56 meV/atom, to be 
compared with 8.3 meV/atom obtained for the fee phase 
at 25 K and low pressure (P = 0.3 GPa, giving a molar 
volume of 9.95 cm'^/mol). This value of Ek is similar 
to those obtained earlier from PIMC simulations in the 
NVT ensemble at temperatures close to 25 K and molar 
volumes around 10 cm'^/mol p3 |. 

At low temperature the ratio E^/Ef. between the ki- 
netic energy of '^He and ^He is close to 1.155, the expected 
value in a harmonic model of lattice vibrations. This in- 
dicates that, irrespective of the important anharmonicity 
present in these solids, anharmonic shifts in Ek scale with 
mass approximately as in a quasiharmonic model. The 
ratio El/E^ decreases as temperature rises and quan- 
tum effects are less important. At the highest simulated 
temperatures (300 K for P = 12 GPa and 600 K for 52 
GPa) we find E^/E^ w 1.05, still clearly larger than the 
classical limit: Ef./E^ — > 1. 
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FIG. 4: Temperature dependence of the kinetic energy of 
hep sohd hehum at two different pressures, as derived from 
PIMC simulations: P = 12 GPa (squares) and 52 GPa (cir- 
cles). Open and filled symbols correspond to *He and '^He, 
respectively. Error bars of the simulation results are less than 
the symbol size. Dashed lines are guides to the eye. 
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FIG. 5: Kinetic-to-potential energy ratio, Ek/Ep, for solid 
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open symbols correspond to fee and hep phases, respectively. 
Error bars, where not shown, are smaller than the symbol 
size. 



IV. DISCUSSION 



The anharmonicity of the lattice vibrations is clearly 
noticeable when one compares the kinetic and potential 
energy of solid helium. From our PIMC simulations, we 
have found in all cases considered here that the vibra- 
tional potential energy Ep is smaller than the kinetic en- 
ergy Ek- Shown in Fig. 5 is the ratio Ek/Ep for "^He 
as a function of pressure at 25 K. We have included in 
this figure results for fee (low pressure) and hep helium 
(high pressure). The kinetic-to-potential energy ratio de- 
creases with increasing pressure; first it gows down fast 
for pressures lower than 5 GPa, and then it decreases 
slower at higher pressures, approaching the value E^/Ep 
= 1, characteristic of harmonic vibrations. We note that 
the fact Ek > Ep has been also found for solids of heavier 
rare gases from PIMC simulations 0, |3^ . 

For solid ^He at low pressure, we find that Ek is larger 
than Ep by about 40%. This important difference be- 
tween Ek and Ep reflects the large anharmonicity of lat- 
tice vibrations in this solid. For comparison, we mention 
that these energies differ by about 20% and 7% for solid 
neon and argon, respectively These numbers 

are very large when compared with covalent or metallic 
solids, and are due to the weakness of van der Waals-type 
bonds, which allows for a large amplitude of the atomic 
vibrations. 



Path-integral Monte Carlo simulations are well suited 
to study quantum effects on structural and thermody- 
namic properties of solids. These effects are particu- 
larly important for solid helium, where isotopic effects 
are relevant, as manifested in differences between crystal 
volume and vibrational energies of solid '^He and ^He. 
The PIMC method enables us to study phonon-related 
properties without the assumptions of quasiharmonic or 
self-consistent phonon approximations, and to study an- 
harmonic effects in solids in a nonperturbative way. This 
method allows us to separate the kinetic and potential 
contributions to the vibrational energy of a given solid, 
and quantify the anharmonicity of the lattice vibrations, 
which together with zero-point motion give rise to iso- 
topic effects. 

PIMC simulations yield in principle "exact" values for 
measurable properties of many-body quantum problems, 
with an accuracy limited by the imaginary-time step 
(Trotter number) and the statistical error of the Monte 
Carlo sampling. Thus, the effectve interatomic potential 
employed here gives a good description of structural and 
thermodynamic properties of solid helium up to pressures 
on the order of those presently achieved in experiments. 
In particular, the equation of state P-V is well described 
by this potential in the pressure range considered here 
(P < 50 GPa), even though it is partially an empiri- 
cal potential, in the sense that the attractive three-body 
part was rescaled ad hoc from the fit to ab initio cal- 
culations. Explicit consideration of three-body terms in 
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PIMC simulations makes the computation times much 
longer, in comparison with simulations including only 
two-body terms. But this has allowed us to check the 
validity limit of some approximations made earlier in this 
kind of calculations, such as including the effect of three- 
body terms in a perturbative way (from configurations 
obtained in PIMC simulations) |ig.ll7l |. As a result, this 
approximation yields at P < 30 GPa an EOS indistin- 
guishable within error bars from that obtained here. At 
higher pressures, however, it predicts volumes larger than 
those found in the present work (and in experiment). 

Due to the large anharmonicity of lattice vibrations in 
solid helium, the isotopic effect on the crystal volume is 
important. At low temperatures, the dependence of the 
crystal volume on isotopic mass is related to the zero- 
point lattice expansion. The magnitude of this isotopic 
effect decreasees appreciably as temperature or pressure 
are raised. Nevertheless, we emphasize that it is still 
measurable at the highest pressures considered here and 
at T close to the corresponding melting temperature, as 
shown in Fig. 3. This is in line with the observed isotopic 
effect on the melting temperature at low pressure [TslllGj . 

Anharmonic effects in solids have been quantified ear- 
lier by comparing the kinetic and potential energies de- 
rived from PIMC simulations [T^ Is^ . For solid helium, 
and for different pressures and temperatures, we have 
found Ek > Ep, as for solids of heavier rare gases with 
Lennard- Jones- type potentials [s^ . The departure of the 
ratio Ek/Ep from its value for a purely harmonic solid 
(Ek/Ep = 1) is a measure of the overall anharmonicity 
of the considered solids. For given T, we have found that 
this energy ratio decreases as pressure is raised. Note 
that the ratio Ek/Ep alone does not give any informa- 
tion on the details of the anharmonicity of the vibra- 
tional modes, but gives a quantitative measure of the 
overall anharmonicity of a given solid. In particular, it 
is well suited to compare the anharmonicity of similar 
materials. For rare-gas solids at low T, Ek/Ep increases 
from 1.02 for xenon to 1.23 for neon at zero pressure [3^ . 
This ratio is clearly higher for solid helium at low P, as 
shown in Fig. 5: Ek/Ep « 1.4 for fee ''He at 25 K. 

It has been recently su gges ted that pressure causes a 
decrease in anharmonicity [S^, l39l l40l | , in line with earlier 
observations that the accuracy of quasiharmonic approxi- 



mations increases as pressure rises and the density of the 
material increases |0, Thus, such approximations 
become exact in the high-density limit js^, • This is 
related to the fact that the ratio of the vibrational energy 
to the whole internal energy of the crystal decreases for 
increasing pressure, in spite of the increase in zero-point 
energy caused by the rise in vibrational frequencies. It 
has been also argued that at high pressures, thermody- 
namic properties of solids can be well described by clas- 
sical calculations, i.e., dealing with the atoms as classical 
oscillators in a given potential This means that, at 
a given temperature, such a classical approach becomes 
more and more accurate as pressure is raised. The origin 
of this is similar to that described above for the success 
of quasiharmonic approaches. Since in this respect the 
lattice vibrations become less relevant as pressure rises, 
and eventually give a negligible contribution to the free 
energy of the solid, their description by a classical or a 
quantum model becomes unimportant for solids in the 
limit of very large pressures. 

In summary, we have carried out PIMC simulations of 
solid helium in the isothermal-isobaric ensemble. In gen- 
eral, our results support the assumptions made in earlier 
calculations, where three-body terms were included in a 
perturbative way. In particular, the equation of state 
obtained here coincides at P < 30 GPa with that found 
in calculations employing the NVT ensemble with that 
assumption. At high pressures, however, data yielded by 
the present PIMC simulations differ from those obtained 
earlier, using similar interatomic potentials. The kind 
of effective interatomic potentials employed here, includ- 
ing two- and three-body terms will probably give a poor 
description of solid helium at very high pressures (> 60 
GPa), but seems appropriate for the pressure range in 
which experimental data are available at present. 
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